Dynamical structure factor of random antiferromagnetic Heisenberg spin chains 
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Combining quantum Monte Carlo simulations with the maximum entropy method, we study the 
dynamical structure factor S(k,uj) of spin- 1/2 antiferromagnetic Heisenberg chains with various 
random bond distributions. We emphasize the crossover behavior in the dynamical properties 
from pure chain to disorder-dominated random singlet phase due to bond randomness. For the 
distribution corresponding to the infinite randomness fixed point, S(k, uj) develops broad non-spinon 
excitations as well as the random-singlet peak near (k, u>) = (n, 0), consistent with the known results 
obtained by the real-space renormalization group method. For weak disorder, however, we find clear 
signature of spinon excitations, reminiscent of pure spin chains, blurred by disorder. We discuss 
the implication for experiments on random-bond antiferromagnetic spin chains, realizable, e.g., in 
BaCu 2 (Si a; Gei-^aCv. 

PACS numbers: 75.10.Nr,75.40.Gb,75.40.Mg 
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I. INTRODUCTION 

The development of the idea of renormalization group 
dramatically advanced our understanding of impurity 
problems^ One branch of the development is the real- 
space renormalization group (RSRG) method,—^ or 
sometimes the strong disorder renormalization group 
method. 5 Originally,^ the RSRG method appeared as 
an algorithm guiding us to approach lower and lower en- 
ergy scales in a controlled fashion, which can be easily 
implemented numerically. Later in a beautiful paper,— 
Fisher demonstrated that it can be firmly laid on analytic 
ground. The central idea of the RSRG method (e.g., in 
a spin system) is at each step we find out the two spins 
coupled with the strongest coupling, replace them with 
an inert spin singlet or a composite spin depending on the 
nature of their coupling, and renormalize the couplings 
in the rest of the system perturbatively. This leads to a 
collection of singlets or other composite spins (seemingly) 
randomly formed at longer and longer length scales. The 
justification of this sequential perturbative approach in a 
random system is that we have or will generate a broad 
energy distribution to ensure its validity. In the three- 
dimensional doped semiconductor casej2£ we have a god- 
send broad energy distribution, because in the insulating 
phase the exchange coupling between local moments (due 
to shallow impurities) decays exponentially with distance 
arising from their hydrogenic envelope functions. We can 
thus begin with the broad distribution, which remains 
broad as the RSRG procedure is carried out. 

However, the situation is not so straightforward in 
disordered spin chains. Quite often, after introducing 
impurities by chemical substitution, we are left with 
a binary distribution of exchange couplings, which can 
be either ferromagnetic or antiferromagnetic in general. 
The neighboring exchange couplings being likely equal 



in strength makes the direct application of the RSRG 
method questionable. This of course is not a matter 
of challenging the low-temperature RSRG behavior of 
such a system, but a step toward understanding physics 
at higher energy scales nevertheless more relevant to 
the interpretation of experimental observations. In the 
case of ferromagnetic chains with diluted antiferromag- 
netic impurity bonds, one of us and collaborators ap- 
plied the modified spin-wave approach to each ferro- 
magnetic segments connected by antiferromagnetic cou- 
plings to capture the higher-energy physics Z The in- 
clusion of these higher-energy excitations beyond "zero 
modes" (which leads to the RSRG results) turns out to 
be important in the description of the crossover behavior 
between high-temperature paramagnetic phase and low- 
temperature random- walk-style, large-spin cluster forma- 
tion predicted by the RSRG analysis. The method may 
not be directly applicable to other problems in general, 
but the idea is general enough and the message is clear: 
we ought to be really careful when deciding which energy 
excitations we need to keep in the renormalization pro- 
cedure (if it is still applicable), especially when we are 
interested in the thermodynamic properties away from 
the so-called infinite disorder fixed point at experimen- 
tally accessible temperature ranges. 

Let us consider the effects of bond randomness in an 
antiferromagnetic (AF) Heisenberg spin chain 
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where Jj > are independent random variables with an 
identical distribution P(J). This model, since studied 
by Dasgupta and Ma£, has been one of the most stud- 
ied examples exhibiting nontrivial effects of quenched 
disorder in a quantum many-body system. Like pure 
spin chains, random spin chains also show spin quan- 
tum number dependent properties. It is well known that 
the ground state of a spin-1/2 AF Heisenberg chain with 
bond randomness flows into the random singlet (RS) 
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phase, which is controlled by an infinite disorder fixed 
point (IRFP) with P*(J) = aJ~ 1+1/Q , no matter how 
weak the strength of disorder is.— ^ For spin-1 AF chain, 
the ground state evolves from the Haldane gapped phase 
at small disorder to the RS phase at strong disorder, with 
a quantum Griffiths phase separating the two.— The cases 
of higher spin quantum number seem to be even more 
complicated as more than one random singlet phase can 
appear i^ilfi Alternative numerical methods, such as quan- 
tum Monte Carlo (QMC) simulatio n 12 i 18 and density ma- 
trix renormalization group (DMRG) method 1 ^ have also 
been implemented in these systems. Besides theoretical 
works, there have been several experimental attempts to 
study the effects of bond randomness in AF Heisenberg 
chainsffi^^ 

Recently, dynamical and transport properties of AF 
spin-1/2 and spin-1 chains have been studied by the 
RSRG method^ In particular, the dynamical structure 
factor for the IRFP is found to be 

S* a P{k,u) = £ $(j q |V2 ln(^H 1/2 ), (2) 
win (il/uj) s 

where a[3 = H — or zz. A and are non-universal 
factors, while $(s) is a universal function. f2 is 
a cut-off in energy scale. This dynamical scaling 
property has been tested by inelastic neutron scatter- 
ing in BaCu2(Sio.5Geo. 5)207 compound. 24 In general, 
BaCu2(Sii_ x Ge x )207 can be regarded as solid solutions 
of BaCu2Si207 and BaCu2Ge207, both prototypical S = 
1/2 Heisenberg AF chains with interchain coupling con- 
stant J = 24.1 meV and 50 meV, respectively. The al- 
ternation of Si and Ge changes the bonding angles of 
neighboring Cu-0 bonds in the zigzag spin chains, due 
to the difference of the covalent radii of Si and Ge. For 
such an almost ideal realization of random spin chain 
with binary exchange coupling distribution, the experi- 
mental group first reported some discrepancies from the 
theoretical prediction for infinite disorder jii One of the 
discrepancies is that the shallow minima surrounding the 
central peak due to dominant AF correlations are miss- 
ing in the experimental plots. The authors suggested 
that the experimental sample may not model an infinitely 
disordered system wellj^ In particular, the higher-order 
renormalizations of the spin operators not captured in 
the theory may affect the shape of the structure factor 
in such a system with finite disorder. The discrepancy 
prompts an interesting question: How does the strength 
of disorder affect the dynamical structure factor? Later, 
through more careful analysis of the original data as 
well as the new data observed for larger samples over 
a larger energy range, the authors of Ref. [2J claimed 
that the momentum-integrated intensity of the dynam- 
ical structure factor measured in BaCu2(Sio. 5660.5)207 
agrees with the result expected by the Miiller ansatz 31 
for a disorder-free spin chain with an assumed coupling 
strength of J — 37 meV. This apparently surprising ob- 
servation deserves more careful investigations of the dy- 
namical properties of disordered spin chains. 



In the experimental compound, the distribution of 
bonds is binary, choosing randomly but with a fixed prob- 
ability between two different coupling strengths, which 
in reality are often not too different from each other in 
value. Apparently, the binary bond distributions P(J) 
is considerably narrower than a power-law distribution 
P*(J) (corresponding to the IRFP). In literature, the 
emergence of such an IRFP distribution from a more re- 
alistic binary distribution is often neglected, because it 
is irrelevant to the fixed-point properties. But in fact, 
the so-called decimation procedure of the RSRG trans- 
formation can fail for the binary distribution initially. 
Though, it is not difficult to amend. In the binary 
case, one can start from AF spin-1/2 segments bonded 
by stronger couplings; these segments are separated by 
weaker bonds. Each segment can thus be renormalized 
by a single spin-1/2 or spin object, depending on the 
number of spins in the segment being odd or even, in 
the same spirit as the normal RSRG procedure. The 
renormalized spin-1/2 objects are effectively coupled, by 
much weaker AF bonds, to one another or to the orig- 
inal spin-1/2 ones surrounded by weak couplings. The 
strength of such a generated bond depends on the length 
of the corresponding spin segments. Obviously, excita- 
tions within the spin segments are hence neglected in 
the subsequent RSRG approach, in many ways similar 
to the case dominated by ferromagnetic bonds X As a 
result, the RSRG approach properly describes the low- 
energy behavior around (k, u) = (w, 0), while leaving the 
behavior at high energies and away from it questionable, 
as addressed by experiments. It is thus desirable to in- 
vestigate the dynamical structure factor to greater detail 
by alternative numerical methods, especially for binary 
distributions. 



In this article, we study the dynamical structure fac- 
tor S(k,w) by combining a QMC approac h 25 i 26 ' 27 with 
the the maximum entropy (Max-Ent) method=&^ in the 
full (k,ui) space. We also use exact diagonalization (ED) 
method in small systems to help us understand the be- 
havior of S(k, u) in random spin chains. We find that for 
the IRFP bond distribution P*(J), the corresponding dy- 
namical structure factor S*(k,u>) is consistent with the 
known theoretical results^ However, like in the above 
mentioned experiment^ we find no evidence of the shal- 
low minima surrounding the central AF correlation dom- 
inated peak predicted by the theory. For weaker disorder, 
we clearly observe a blurred lower bound that follows the 
spinon dispersion relation, in addition to the weakened 
low-energy behavior of S(k,uj) around (fc, u>) = (jr,0). 
As expected, we find the strength of the random bond 
distribution affects the properties of S(k,u) in the (fc, to) 
space away from (k,ui) = (tt, 0), where the dynamical 
behavior is dominated by higher-energy excitations not 
included in the RSRG approach. After we present our 
main results in Sec. HU we will discuss the implication of 
our results to experiments in Sec. IIIII 



II. NUMERICAL RESULTS 

The QMC simulation with continuous imaginary time 
loop/cluster algorithm 26 has been successfully applied to 
bond disordered system^ 2 . On the other hand, the Max- 
Ent method, which analytically continues Euclidean-time 
QMC data to real frequencies, has been successfully ap- 
plied to spin- 1/2 and spin-1 AF Heisenberg chains at low 
and high temperature o 13 ' 14 i 28 with or without random- 
ness. In this work, we adapt the Max-Ent implemen- 
tation by Jarrell and Gubernatis.^ We have tested our 
QMC simulation with the Max-Ent method on various 
models, include spin-1/2 (spin-1) AF Heisenberg chains 
and spin-1/2 AF chain with a magnetic field, and ob- 
tained consistent results with those in Ref. [28|, as well 
as in Ref. UtI For random antiferromagnetic Heisenberg 
chains, we simulate systems with up to L — 128 spins and 
inverse temperature (3 = 100, in units of the averaged 
bond strength J. We study several random-bond distri- 
butions P(J), including (1) the IRFP distribution P*(J) 
with a = l.OiS, (2) the flat distribution of J e [1,2], 
which we label as Pf and (3) the binary distributions, 
where the bond strength is chosen to be J = 2.0 with 
probabilities P b = 0.1 — 0.5 and J = 1.0 with probabil- 
ity (1 — P&). For each P(J), we average over 800 real- 
izations, and simulate each realization with 2000 Monte 
Carlo sweeps (MCS) for thermalization due to the dy- 
namical anisotropy in disordered system^, 400 MCS for 
measurements to be binned for the Max-Ent method to 
remove correlations of measurements. We calculate the 
imaginary-time correlation function by 
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G(k,r) = Y,Y, eik{i ~ l) < sf(r)s;(o) > 



(3) 
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with 100 slices in imaginary time. S(k,u>) are then ex- 
tracted by the Max-Ent method. The application to the 
bond random quantum lattice model is straightforward. 
First, for each sample of bond random realizations i, 
Max-Ent method can give its corresponding dynamical 
spectra Si(k,w), and the statistical average < Si(k,co) > 
of bond random realizations can be easily obtained. Sec- 
ond, as each bond random realization is statistically in- 
dependent, disorder effect can be extracted because the 
Bayesian approach to data analysis incorporates prior 
knowledge about probabilities relationships among the 
data. 

We also calculate magnetic susceptibility Xu{T) as a 
function of temperature T. Figure [T] shows that \ u for 
two different P(J) approaches Curie like behavior at low 
temperatures. The feature is consistent with the typical 
behavior of the RS phase 
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We can extract the energy cutoff f2 from the numerical 
fit of Xu{T) to Eq. (g|), and the results are Q = 2.6 for 
P*{J) and fi = 2.2 for P b = 0.5. 
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FIG. 1: Log-log plot of magnetic susceptibility Xu(T) vs. 
temperature T. For comparison, we also plot a Curie law 
Xu oc 1/T (dashed straight line). 
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FIG. 2: S(k,cj) at fixed w = 0.1J of P*(J) with a = 1. The 
filled points are results of the combined QMC and Max-Ent 
method, and the solid line is the best fit to the analytical 
result of RSRG [Eq. ©] near k = tt. 



We first calculate S(k,u>) by the combined QMC and 
Max-Ent methods for the IRFP distribution P*(J) with 
a = 1.0, for which the analytical result [Eq. ([2])] has been 
known in the RSRG approach. Figure [2] shows S(k, oj) 
for uo = 0.1J. The fluctuations of the QMC data are 
due to the intrinsic bias of the Max-Ent method used for 
the analytical continuation. We can fit the data in the 
vicinity of (k— tt) <C 1 to Eq. @, and one can see that our 
results are consistent with theoretical predictions. We 
note that, as pointed by the authors of Ref. HU, we do 
not observe the pair of shallow minima in S(k, to) around 
(fc,w) = (7r,0). 

We show our results of S(k,uj) in 3D plots for the 
IRFP distribution P*(J) (Fig. [3]), the flat distribution P f 
(Fig. g}, the binary distribution with P b = 0.1 (Fig. [5J), 
and the binary distribution with P b = 0.5 (Fig. [S]). The 
low-energy excitation peak at (k,oS) = (tt, 0) is visible 
in all these figures, indicating the antiferromagnetic na- 
ture of the couplings. For the IRFP distribution P*(J), 
there are a considerable amount of low-energy excitations 
below the lower boundary of the two-spinon continuum 
(in the corresponding pure chain) uji = sin(fc) and 
also high-energy excitations above the upper boundary 
of oj u = 7rJsin(-|), where J is the arithmetic average of 
bonds. For the flat distribution Pf, we note one can still 
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FIG. 3: S(jfe,o>) for the IRFP distribution P*(J) for a = 1.0. 
The solid lines are the boundary of the two-spinon continuum 
for a pure spin- 1/2 AF Heisenberg chain with average of bonds 
J = 0.5. The palette on the right gives the strength of S(k, ui). 
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FIG. 4: S(k,uj) for the flat distribution P f for J G [1, 2]. The 
solid lines are the boundary of the two-spinon continuum for 
a pure spin-1/2 AF Heisenberg chain with average of bonds 
J= 1.5. 



recognize the blurred boundaries of the two-spinon con- 
tinuum. This suggests that spinon-like excitations dom- 
inate the excitations in the system and there are only 
limited excitations beyond these boundaries. The results 
reveals that under weak disorder strength, such as for the 
flat distribution, the low-energy excitations can still be 
described by spinons with a broadened dispersion curve 
(due to the weak fluctuation of bond strength) . The re- 
sults of the binary distribution with Pb = 0.1 is very 
similar as those of the flat Pf. In this case, it is easier to 
imagine the elementary excitations are spinons for spin- 
1/2 AF spin segments of various lengths (weakly coupled 
together). 

For the Pb — 0.5 binary distribution, the boundaries of 
the two-spinon continuum seem to be blurred more than 
those for the flat Pf and Pb = 0.1 binary distributions, 
but certainly less than the IRFP distribution P*(J). It 
suggests that the disorder strength for the Pb — 0.5 bi- 
nary distribution can be considered as in between the flat 
Pf or Pb = 0.1 binary distribution and the IRFP distri- 
bution P*(J). It is not difficult to quantify the disorder 
strength thus explain the results by defining the variance 
of bond randomness 



^(ln^-Jh^J (5) 

where J is the arithmetic average of bonds. As expected, 
we obtain <5 = 1.0 for P*(J), S = 0.346 for P b = 0.5, 
6 = 0.208 for P b = 0.1, and 5 = 0.197 for P f . 

To further compare S(k,u>) for different distributions, 
we plot our results at oj = 0.1 J, 0.3 J, and 0.5 J in Fig. [7] 



For oj = 0.1 J, S(k,oj) peaks at (k,ui) = (7r,0) for the 
IRFP distribution P*(J). However, S(k,ui) peaks at 
k < 7T for the weaker distributions, due to the reminiscent 
of spinon excitations. At a higher to — 0.3 J, S(k,uj) still 
peaks at k = it for the IRFP distribution P*(J) but the 
peak is much broader than at lu = O.U. The spinon exci- 
tation peak in S(k, lu) for the flat Pf and Pb = 0.1 binary 
distributions shifts significantly away from k — n, while 
S(k,Lu) is almost flat for the Pb = 0.5 binary distribu- 
tion, reflecting the interplay of the IRFP and the spinon 
excitations in the pure case. At still higher u>, S(k,Lu) 
peaks at values far away from k = tt for all distributions. 

In all our calculations by Max-Ent method for different 
models, we find the qualitative features of S(k, uo) remain 
changed. Because of the strong correlations of different to 
in the Max-Ent method, the specific error bar for each ui 
lacks meaning, and it is only possible to assign error bars 
to the integrated functions of S(k,u)£2- Thus we don't 
plot them in all related figures in this paper. However we 
note that all error bars are within acceptable numerical 
precision in all these figures. 

To relate the strength of disorder to a measurable 
quantity, we introduce the normalized weight of S(k, to) 
confined outside the spinon continuum region 

where T is the spinon continuum region T — {(fc,u;)|fc £ 
[0, 7r],w 6 [n^ sin(fc), 7r J sin(fc/2)]}, and 5 stands for the 
entire region of {(k,ui)\k e [0, tt],u) £ [0, oo]}. Not sur- 



k 



FIG. 5: S(k,u>) for the binary distribution with Pt = 0.1 
(i.e. 90% J = 1 and 10% J = 2). The solid lines are the 
boundary of the two-spinon continuum for a pure spin- 1/2 
AF Heisenberg chain with average of bonds J = 1.1. 

prisingly, we find r roughly increases with 8. We note 
that a quantity based on the similar idea has already 
been used to analyze the disorder effect on the dynami- 
cal structure factor of transverse- field Ising chains. 30 



III. SUMMARY AND DISCUSSION 

In this paper, we combined QMC simulations with the 
Max-Ent method to investigate the dynamics of spin-1/2 
AF Heisenberg chains with various random bond distri- 
butions. Our results are consistent with the previous the- 
oretical predictions, including the IRFP-controlled be- 
havior. While all these systems go into the RS phase 
irrespective of the random strength, the strength of the 
bond random distribution plays an important role in the 
crossover region from a pure AF chain to the RS phase. 
For weak randomness, spinon excitations can dominate 
the dynamical structure factor S(k 1 to) much like in the 
pure case. 

Our results also reveal that for the quasi-one- 
dimensional materials with impurity bonds, the behavior 
of S(k, lo) can look very similar as the pure case. On the 
other hand, non-spinon excitations, visible in the simu- 
lation results, may be buried in the background signal 
in experiments, subtracted when analyzing inelastic neu- 
tron scattering data. For BaCu2(Sio. 5660.5)207, Masuda 
et al£^ noted that even for the relative strong bond ran- 
domness Pb = 0.5, the effects of disorder are difficult to 
extract from inelastic neutron scattering data, somewhat 
different from their earlier claim. To further clarify the 
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FIG. 6: S(k,ui) for the binary distribution with Pb = 0.1 
(i.e. 50% J = 1 and 50% J = 2). The solid lines are the 
boundary of the two-spinon continuum of a pure spin-1/2 AF 
Heisenberg chain with average of bonds J = 1.5. 

issue, we show our results using the same scaling rela- 
tional as these authors did in Fig. [9l According to the 
Muller ansata 3 ! for the pure chain 

s(fc,«) = -^=e(w-w,)eK-w), (7) 

v w - ^ 

ujS(u) — w J™ S(k,ui)dk is linear in ui at low energies. In 
the above equation, 0(x) is the usual step function. We 
found that this remains to be true for Pi, = 0.1, 0.3, and 
0.5. For example for P b — 0.1, ujS(u>) increases roughly 
linearly with ui in low energy region, consistent with the 
Muller ansatz result. We note there are small fluctu- 
ations around a straight line, which is induced by the 
statistical bias of Max-Ent method visible in all these 
ujS{lo) curves in Fig. [5J The results are in good agree- 
ment with the claim in the Erratum in Refs.[2J. We note 
that for Pt = 0.5, wS(w) bends down at w « 0.6J in 
our simulational result, implying the breakdown of the 
Muller ansatz. A recent experimented has reported a 
breakdown at a similar energy scale. This breakdown 
can be attributed, apparently, to disorder, as the curves 
for smaller doping concentrations bend down less at high 
energy scales. 

For a binary distribution of bond strength J\ and J2 
( J\ < J2) with probability (1 — P&) and P&, respectively, 
the disorder strength [defined in Eq. (JSJ)] is found to be 

S 2 = P b (l-P b )\n 2 (J 2 /J 1 ). (8) 

5 reaches the largest value at P b = 0.5. Thus one should 
see even weaker disorder effect for x > 0.5 or x < 0.5 in 
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FIG. 9: ojS(w) for various binary distributions. Black circles 
are data of Pt=0.5, red squares are data of Pt=0.3, and blue 
triangles are data of P b =0.1. 



BaCu 2 (Sii_ 2; Ge 2: )207. Therefore, we conclude the effects 
of quench disorder may only be seen in strong disorder 
case experimentally difficult to realize. Of course, some 
of the experimental difficulties, as pointed out, e.g., in 
Ref. HH, lie on the experimental wave vector resolution, 
which relates to the lowest energy transfer. On the other 
hand, our suggestion of measuring the normalized weight 
of SlkjUj) confined outside the spinon continuum region 
can be an indication, but only qualitative, due to under- 
lying broad background noised It remains interesting 
to see whether the disorder effect can be observed in the 
dynamical structure factor (as measured by the normally 
informative inelastic neutron scattering experiments) for 
a realistic quasi-one-dimensional system. 



FIG. 7: S(k, u) at u = 0.1J, 0.3 J, and 0.5J. The black circles 
are results for IRFP distribution P* ( J), red squares for binary 

distribution Pu =r 5 orppn fillpd diamonds for R. = D 1 and 
blU! 
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FIG. 8: Normalized weight of S(k,u>) outside the spinon con- 
tinuum as function of S. 



Acknowledgments 



We thank Kun Yang for very helpful discussions. This 
work is supported by the National Natural Science Foun- 
dation of China through Projects 90303010 (H.P.Y.) and 
10504028 (X.W.). 



1 K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975). 



2 C. Dasgupta and S.-K. Ma, Phys. Rev. B 22, 1305 (1980). 



7 



3 R. N. Bhatt, and P. A. Lee, Phys. Rev. Lett. 48, 344 
(1982). 

4 D. S. Fisher, Phys. Rev. B 50, 3799 (1994); ibid. 51, 6411 
(1995). 

5 F. Igloi and C. Monthus, Phys. Rep. 412, 277 (2005). 

8 K. Andres, R. N. Bhatt, P. Goalwin, T. M. Rice, and R. E. 
Walstedt, Phys. Rev. B 24, 244 (1981). 

7 X. Wan, K. Yang, and R. N. Bhatt, Phys. Rev. B 66, 
014429 (2002). 

8 R. A. Hyman, and K. Yang, Phys. Rev. Lett. 78, 1783 
(1997). 

9 G. Refael, S. Kehrein, and D. S. Fisher, Phys. Rev. B 66, 
060402 (2002). 

10 K. Damle and D. A. Huse, Phys. Rev. Lett. 89, 277203 
(2002). 

11 K. Damle, O. Motrunich, and D. A. Huse, Phys. Rev. Lett. 
84, 3434 (2000); O. Motrunich, K. Damle, and D. A. Huse, 
Phys. Rev. B 63, 134424 (2001). 

12 S. Todo, K. Kato, and H. Takayama, J. Phys. Soc. Jpn. 
Suppl. A 69, 355 (2000); T. Arakawa, S. Todo, and H. 
Takayama, J. Phys. Soc. Jpn. 74, 1127 (2005). 

13 A. W. Sandvik, Phys. Rev. B 52, R9831 (1995); O. A. 
Starykh, A. W. Sandvik, and R. R. P. Singh, ibid. 55, 
14953 (1997). 

14 S. Wessel and S. Hass, Phys. Rev. B 65,132402 (2002). 

15 The fortran code can be download from 



|http://ww.physics.uc.edu/~Jarrell/LINUX/MEMj 

16 We note, strictly speaking, the IRFP distribution should 
be referred to the limit a — > oo. We use the term here, 
nevertheless, for convenience. 

17 M. B. Stone, D. H. Reich, C. Broholm, K. Lefmann, C. 
Rischel, C. P. Landee, and M. M. Turnbull, Phys. Rev. 
Lett. 91, 037205 (2003). 



N. Laflorencie, H. Rieger, A. W. Sandvik, and P. Henelius, 

Phys. Rev. B 70, 054430 (2004) 

K. Hida, Phys. Rev. Lett. 83, 3297 (1999). 

J. C. Scott, A. F. Garito, A. J. Heeger, P. Nannelli, and 

H. D. Gillman, Phys. Rev. B 12, 356 (1975). 

Y. Endoh, G. Shirane, R. J. Birgeneau, and Y. Ajiro, Phys. 

Rev. B 19, 1476 (1979). 

L. C. Tippie, and W. G. Clark, Phys. Rev. B 23, 5846 
(1981). 

J. C. Boucher, C. Dupas, W. J. Fitzgerald, K. Knorr, and 
J. P. Renard, J. Phys. (Paris) 39, L86 (1978). 
T. Masuda, A. Zheludev, K. Uchinokura, J. H. Chung, 
and S. Park, Phys. Rev. Lett. 93, 077206 (2004); ibid. 96, 
169908 (2006). 

H. G. Evertz, Adv. Phys. 52, 1 (2003). 

U. J. Wiese, and H. P. Ying, Z. Phys. B 93, 147 (1994); 

B. B. Beard, and U. J. Wiese, Phys. Rev. Lett. 77, 5130 

(1996). 

N. Kawashima, and K. Harada, J. Phys. Soc. Jpn. 73, 1379 
(2004). 

R. N. Silver, D. S. Sivia, and J. E. Gubernatis, Phys. Rev. 
B 41, 2380 (1990); J. E. Gubernatis, M. Jarrell, R. N. 
Silver, and D. S. Sivia, ibid. 44, 6011 (1991). 
M. Jarrell and J. E. Gubernatis, Phys. Rep. 269, 133 
(1996); J. Deisz, M. Jarrell, and D. L. Cox, Phys. Rev. 
B 42, 4869 (1990); J. Deisz, M. Jarrell and D. L. Cox, 
ibid. 48, 10277 (1993). 

X. Jia and S. Chakravarty, Phys. Rev. B 74, 172414 (2006). 
G. Miiller, H. Thomas, H. Beck, and J. C. Bonner, Phys. 
Rev. B 24, 1429 (1981). 

A. Zheludev, T. Masuda, G. Dhalenne, A. Revcolevschi, C. 



Frost, and T. Perring, cond-mat/0611030 (unpublished). 



